#########################################################################   
#                                 INFO                                  #   
#########################################################################  

  # PROJECT: Gender quotas in Tunisia's 2018 municipal elections
  # PURPOSE: Analysis for FHL placement as DV
  # CREATED: December 2018 by Julia Clark
  # INPUTS:  quotas_mun_clean.R, quotas_list_clean.R,
  # OUTPUTS: graphs and tables

#########################################################################   
#               SETUP                                                   #   
######################################################################### 

  ######## ENVIRONMENT
  
  rm(list = ls()) 
  setwd("~/Desktop/replication_what_men_want/")

  ######## PACKAGES
  
  need <- c("openxlsx", "tidyr", "dplyr", "ggplot2", "ggthemes", "extrafont", 
            "margins", "multiwayvcov", "lmtest", "stargazer", "viridis", "miceadds",
            "ggResidpanel", "modelsummary") 
  have <- need %in% rownames(installed.packages()) 
  if(any(!have)) install.packages(need[!have]) 
  invisible(lapply(need, library, character.only=T)) 
  
  ######## SET ALL TABLES TO SHOW "NA" BY DEFAULT
  
  table = function (..., useNA = 'ifany') base::table(..., useNA = useNA)  
  
  ######## CLUSTERING FUNCTIONS
  
  f.test <- function(model, cluster){
    vcovCL<-cluster.vcov(model, cluster)
    ftest <- waldtest(model, vcov = vcovCL, test = "F")
    return(ftest)
  } 
  
  cluster.se <- function(model, cluster){
    coeftest(model, vcov=function(x) cluster.vcov(x, cluster))[,2]
  }

#########################################################################   
#               GGPLOT PARAMETERS                                       #   
#########################################################################    
  
  # Font sizes
  bold16 <- element_text(face = "bold", color = "black", size = 16)
  bold18 <- element_text(face = "bold", color = "black", size = 18)
  bold20 <- element_text(face = "bold", color = "black", size = 20)
  greyitalic18 <- element_text(face = "italic", color = "grey40", size = 18)
  
  # Continuous axes
  cy <- scale_y_continuous(expand = c(0,0))
  cx <- scale_x_continuous(expand = c(0,0)) 
  
  # Colors 
  vid_low <- "#481567FF"
  vid_med <- "#29AF7FFF"
  vid_high <- "#FDE725FF"
  
  # General theme
  lecs_theme <- theme_bw() + 
    theme(text = element_text(size = 20),
          title = bold18, 
          plot.subtitle = greyitalic18,
          axis.title = bold16, 
          legend.title = bold16,
          plot.caption = element_text(face = "plain", color = "grey40", size = 14,
                                      hjust = 0, margin = margin(t = 15, r = 0, b = 0,
                                                                 l = 0)),
          axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)),
          axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0),
                                      angle = 0, vjust = .5),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank()
    )
  
  # Theme with lines for coeff plots
  lecs_theme_lines <- theme_bw() + 
    theme(text = element_text(size = 20),
          title = bold18, 
          plot.subtitle = greyitalic18,
          axis.title = bold16, 
          legend.title = bold16,
          plot.caption = element_text(face = "plain", color = "grey40", size = 1,
                                      hjust = 0, margin = margin(t = 15, r = 0, b = 0,
                                                                 l = 0)),
          axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0)),
          axis.title.y = element_text(margin = margin(t = 0, r = 15, b = 0, l = 0),
                                      angle = 0, vjust = .5),
          panel.grid.minor = element_blank(),
          panel.background = element_blank()
    )
  
  # Smaller margins
  zero_margins <- theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) 

#########################################################################   
#               FUNCTIONS FOR GRAPHING                                  #   
#########################################################################   
  
  ##### PREDICTED VALUES GRAPHS
  
  pred.fun <- function(data, strength, mod, name, turnout, pop, unemploy, urban){
    
    # Set percentiles for variables [remove most extreme 5%]
    min <- round(quantile(data[[strength]], 0.025, na.rm = T), 0)
    max <- round(quantile(data[[strength]], 0.975, na.rm = T), 0)
    range <- seq(min, max, .01)
    
    # Create dataframes
    df <- data.frame(list_class_ennahdha = c(rep("Ennahdha", length(range)*3),
                                             rep("Nidaa", length(range)*3),
                                             rep("Third Party", length(range)*3)),
                     strength = rep(range, 9),
                     strengthsq = rep(range^2, 9),
                     log_pop = pop, 
                     emp_unemploy_rate = unemploy,
                     pop_com_per = urban, 
                     turnout_per_pop_f_2014_parl = turnout,
                     gov_en = "Ben Arous", # Medenine, Siliana, El Kef, 
                     model = name
    ) 
    df <- df %>% filter(list_class_ennahdha %in% data$list_class_ennahdha) 
    
    # Rename strength measure column
    names(df)[names(df) == "strength"] <- strength
    names(df)[names(df) == "strengthsq"] <- paste0(strength,"sq")
    
    # Predict values
    prediction <- predict(mod, df, type = "response", se.fit = T)
    
    # Bind together
    df1 <- cbind(var = names(df)[2],
                 val = range,
                 pred = prediction$fit,
                 se = prediction$se.fit,
                 pred_l = prediction$fit - prediction$se.fit*1.96,
                 pred_h = prediction$fit + prediction$se.fit*1.96,
                 df)
  }
  
  ###### COEFFICIENT PLOTS
  
  # Set variable names
  cm <- c(
    'votes_per_nat_mean_2014_parl' = 'Relative voteshare',
    'votes_per_nat_mean_2014_parlsq' = 'Relative votershare sq',
    'log_pop' = 'Log population',
    'pop_com_per' = 'Urbanization Rate',
    'emp_unemploy_rate' = 'Unemployment Rate',
    'turnout_per_pop_f_2014_parl' = '2014 female turnout',
    'list_class_ennahdhaNidaa' = 'Nidaa',
    'list_class_ennahdhaThird Party' = 'Third Party',
    'list_class_ennahdhaEnnahdha' = 'Ennahdha',
    'votes_per_nat_mean_2014_parl:list_class_ennahdhaNidaa' = '2014 votershare centered (nat) x Nidaa',
    'list_class_ennahdhaNidaa:votes_per_nat_mean_2014_parlsq' = '2014 votershare centered sq (nat) x Nidaa',
    'votes_per_nat_mean_2014_parl:list_class_ennahdhaThird Party' = '2014 votershare centered (nat) x Third Party',
    'list_class_ennahdhaThird Party:votes_per_nat_mean_2014_parlsq' = '2014 votershare centered sq (nat) x Third Party'
  )
  
  
#########################################################################   
#               LOAD DATA                                               #   
#########################################################################  

  # Municipal and list-level formatted data
  load("data/quotas_mun_clean.RData")
  load("data/quotas_list_clean.RData")

#########################################################################   
#               SUBSET DATA FOR ANALYSIS                                #   
#########################################################################   
  
  ##### MAIN MODELS
  
  # Complete cases for core variables  
  data.na <- list %>%
    filter(  !is.na(votes_per_nat_mean_2014_parl) & !is.na(votes_per_nat_mean_2014_parlsq) &
             !is.na(turnout_per_pop_f_2014_parl) & !is.na(log_pop) & !is.na(pop_com_per) &
               !is.na(emp_unemploy_rate))
  
  nrow(data.na) # 988
  
  # Use standardized version of variables
  data.std <- data.na %>%
    mutate(votes_per_nat_mean_2014_parl = sd_votes_per_nat_mean_2014_parl ,
           votes_per_nat_mean_2014_parlsq = sd_votes_per_nat_mean_2014_parlsq,
           votes_per_2014_parl = sd_votes_per_2014_parl,
           votes_per_2014_parlsq = sd_votes_per_2014_parlsq,
           votes_per_dis_mean_2014_parl = sd_votes_per_dis_mean_2014_parl ,
           votes_per_dis_mean_2014_parlsq + sd_votes_per_dis_mean_2014_parlsq ,
           list_margin3_2014_parl = sd_list_margin3_2014_parl ,
           list_margin3_2014_parlsq = sd_list_margin3_2014_parlsq ,
           turnout_per_pop_f_2014_parl = sd_turnout_per_pop_f_2014_parl ,
           log_pop = sd_log_pop ,
           emp_unemploy_rate = sd_emp_unemploy_rate ,
           di_relig_pol_index = sd_di_relig_pol_index ,
           di_women_pol_index = sd_di_women_pol_index ,
           sd_fem_per = sd_sd_fem_per ,
           edu_f_high = sd_edu_f_high ,
           edu_m_high = sd_edu_m_high,
           pop_com_per = sd_pop_com_per,
    )
  
  # Set data for main models
  data.main <- data.std
  nrow(data.main) # 988
  
  ##### ROBUSTNESS CHECKS
  
  # Ennahdha and Nidaa only 
  data.nn <- data.main %>% filter(list_class == "Ennahdha" | list_class == "Nidaa") %>%
    mutate(list_class_ennahdha = factor(list_class_ennahdha))
  
  # Biggest strongholds [Tataouine = Ennadhdha, Monastir = Nidaa]
  strong <- data.main %>%
    group_by(gov_en, list_class_ennahdha) %>%
    summarise(mean = mean(votes_per_nat_mean_2014_parl, na.rm = T))
  data.nostrong <- filter(data.nn, gov_en != "Tataouine" & gov_en != "Monastir")
  
  # Dropping largest cities
  data.nolarge <- data.main %>%
    filter(mun_en != "Tunis" & mun_en != "Sfax" & mun_en != "Sousse" )
  nrow(data.nolarge) # 977
  
  ##### PARTY-LEVEL ANALYSIS
  
  # Individual classes
  data.en <- data.main %>% filter(list_class == "Ennahdha")
  data.ni <- data.main %>% filter(list_class == "Nidaa")
  data.tp <- data.main %>% filter(list_class == "Third Party")
  
#########################################################################   
#               FHL PLACEMENT: MAIN RESULTS   (TABLE 1)                 #   
#########################################################################      
  
  # THEORY: Parties will place FHLs strategically in municipalities 
  # that will help them preserve important incumbents at the heads of 
  # lists. This implies that FHLs will be *most* likely where: 
  
    # H1. The party is relatively weak [2014 relative voteshare]
    # H2. The importance of winning is lower [smaller municipalities]

  ##### OLS MODELS [FOR TABLE 1] #####
  
  # Create list of model names
  strength.ols <- list()
  
  # Run models
  strength.ols[[1]] <- lm(FHL ~ votes_per_nat_mean_2014_parl  + 
                            log_pop + gov_en, data = data.main)
  strength.ols[[2]] <- lm(FHL ~ votes_per_nat_mean_2014_parl  + 
                            log_pop + pop_com_per + emp_unemploy_rate + 
                            turnout_per_pop_f_2014_parl + 
                            list_class_ennahdha + gov_en, data = data.main)
  strength.ols[[3]] <- lm(FHL ~ votes_per_nat_mean_2014_parl + votes_per_nat_mean_2014_parlsq + 
                            log_pop + pop_com_per + emp_unemploy_rate + 
                            turnout_per_pop_f_2014_parl + 
                            list_class_ennahdha + gov_en, data = data.main)
  strength.ols[[4]] <- lm(FHL ~ above_mean_2014_parl + 
                            log_pop + pop_com_per + emp_unemploy_rate + 
                            turnout_per_pop_f_2014_parl + 
                            list_class_ennahdha + gov_en, data = data.main)
  names(strength.ols) <- c("No controls","Linear", "Quadratic", "Binary")
  
  # Cluster standard errors by municipality
  strength.ols.cl <- list()
  for(i in 1:length(strength.ols)){
    strength.ols.cl[[i]] <- cluster.se(strength.ols[[i]], data.main$u_id)
  }
  
  # Calculate new F-stat with clustering
  check.f <- list()
  strength.ols.f <- NA
  strength.ols.pval <- NA
  for(i in 1:length(strength.ols)){
    check.f[[i]] <- f.test(strength.ols[[i]], data.main$u_id)
    strength.ols.f[i] <- round(check.f[[i]]$F[2],3)
    strength.ols.pval[i] <- round(check.f[[i]]$`Pr(>F)`[2],4)
  }
  print(strength.ols.pval)

  # Output table
  stargazer(strength.ols, 
             out = "tables/fhl_placement_ols.tex",
            title = "Strategic Placement of Female-Headed-Lists (FHLs)",
            label = "fhl.ols",
            table.placement = "H",
            type = "text", 
            column.labels = names(strength.ols),
            dep.var.labels = "Probability of an FHL",
            omit = c("gov_en", "list_class"),
            covariate.labels = c( "Relative strength in 2014",
                                  "Relative strength in 2014 squared",
                                  "Above mean relative strength in 2014",
                                  "Log population",
                                  "Urbanization rate",
                                  "Unemployment rate",
                                  "Female turnout in 2014",
                                  "Constant"),
            se = strength.ols.cl, 
            omit.stat = "f", 
            font.size = "scriptsize",
            add.lines = list(c("Party Controls", "N","Y", "Y", "Y"),
                             c("Governorate FE", "Y","Y", "Y", "Y"),
                             c("F Statistic", paste(strength.ols.f, "***",sep = ""))),
            notes.align = "l",
            notes = c("OLS models with standardized continuous variables (2014 relative strength, log population, ",
                      "urbanization rate, unemployment rate, and female turnout) and controls for party. Standard ",
                      "errors clustered by municipality. Relative strength is the party's voteshare in the 2014 ", 
                      "parliamentary elections at the municipal level minus the party's average voteshare across all", 
                      "municipalities. Analysis excludes independents lists and party lists that ran in 2018 but not 2014.")
      ) 



  ##### LOGIT MODELS [FOR GRAPHING + APPENDIX TABLE A.7] #####
  
  # Create list of model names
  strength.log <- list()
  
  # Run models
  strength.log[[1]] <- glm(FHL ~ votes_per_nat_mean_2014_parl  + 
                            log_pop + gov_en, data = data.main, family = "binomial")
  strength.log[[2]] <- glm(FHL ~ votes_per_nat_mean_2014_parl  + 
                             log_pop + pop_com_per + emp_unemploy_rate + 
                             turnout_per_pop_f_2014_parl + 
                             list_class_ennahdha + gov_en, data = data.main, family = "binomial")
  strength.log[[3]] <- glm(FHL ~ votes_per_nat_mean_2014_parl + votes_per_nat_mean_2014_parlsq + 
                             log_pop + pop_com_per + emp_unemploy_rate + 
                             turnout_per_pop_f_2014_parl + 
                             list_class_ennahdha + gov_en, data = data.main, family = "binomial")
  strength.log[[4]] <- glm(FHL ~ above_mean_2014_parl + 
                             log_pop + pop_com_per + emp_unemploy_rate + 
                             turnout_per_pop_f_2014_parl + 
                             list_class_ennahdha + gov_en, data = data.main, family = "binomial")
  
  # Cluster standard errors by municipality
  strength.log.cl <- list()
  for(i in 1:length(strength.log)){
    strength.log.cl[[i]] <- cluster.se(strength.log[[i]], data.main$u_id)
  }
  
  ###### FHL PLACEMENT: MAIN MODELS ######
  
  # Logit table (with full variables)
  stargazer(strength.log, 
            out = "tables/fhl_placement_log.tex",
            title = "Logit Models: Strategic Placement of Female-Headed-Lists (FHLs)",
            label = "fhl.log",
            table.placement = "H",
            type = "text", 
            column.labels = names(strength.log),
            dep.var.labels = "Probability of an FHL",
            covariate.labels = c( "Relative strength in 2014",
                                  "Relative strength in 2014 squared",
                                  "Above mean relative strength in 2014",
                                  "Log population",
                                  "Urbanization rate",
                                  "Unemployment rate",
                                  "Female turnout in 2014",
                                  "Constant"),
            se = strength.log.cl, 
            omit = c("gov_en","list_class"),
            font.size = "scriptsize",
            add.lines = list(c("Governorate FE","N", "Y", "Y", "Y")),
            notes.align = "l",
            notes = c("OLS models with standardized continuous variables (2014 relative strength, log population, ",
                      "urbanization rate, unemployment rate, and female turnout) and controls for party. Standard ",
                      "errors clustered by municipality. Relative strength is the party's voteshare in the 2014 ", 
                      "parliamentary elections at the municipal level minus the party's average voteshare across all", 
                      "municipalities. Analysis excludes independents lists and party lists that ran in 2018 but not 2014.")
  ) 
  
  
  
  ##### PREDICTED VALUES GRAPHS #####
  
  # Main Models
  pred.main.log <- pred.fun(data = data.main, strength = "votes_per_nat_mean_2014_parl", 
                            turn = mean(data.main$turnout_per_pop_f_2014_parl, na.rm = T),
                            pop = mean(data.main$log_pop, na.rm = T),
                            urban = mean(data.main$pop_com_per, na.rm = T),
                            unemploy = mean(data.main$emp_unemploy_rate, na.rm = T),
                            mod = strength.log[[2]], name = "linear") 

  
  # Population Models
  pred.pop.low <- pred.fun(data = data.main, strength = "votes_per_nat_mean_2014_parl", 
                           turn = mean(data.main$turnout_per_pop_f_2014_parl, na.rm = T),
                           urban = mean(data.main$pop_com_per, na.rm = T),
                           unemploy = mean(data.main$emp_unemploy_rate, na.rm = T),
                           pop = -1,
                           mod = strength.log[[2]], name = "Low (-1 SDs)") 
  
  pred.pop.med <- pred.fun(data = data.main, strength = "votes_per_nat_mean_2014_parl", 
                           turn = mean(data.main$turnout_per_pop_f_2014_parl, na.rm = T),
                           urban = mean(data.main$pop_com_per, na.rm = T),
                           unemploy = mean(data.main$emp_unemploy_rate, na.rm = T),
                            pop = 0,
                            mod = strength.log[[2]], name = "Medium (0 SDs)") 
  
  pred.pop.hig <- pred.fun(data = data.main, strength = "votes_per_nat_mean_2014_parl", 
                           turn = mean(data.main$turnout_per_pop_f_2014_parl, na.rm = T),
                           urban = mean(data.main$pop_com_per, na.rm = T),
                           unemploy = mean(data.main$emp_unemploy_rate, na.rm = T),
                            pop = 1,
                            mod = strength.log[[2]], name = "High (+1 SDs)") 
  
  pred.pop <- rbind(pred.pop.low, pred.pop.med, pred.pop.hig)
  
  # Change types to factor
  pred.pop <- mutate(pred.pop, model = factor(model, levels = c("High (+1 SDs)", "Medium (0 SDs)", "Low (-1 SDs)")))
  
  # Subset data for graphing
  graph.linear <- filter(pred.main.log, list_class_ennahdha == "Ennahdha")
  # graph.linear.lines <- filter(graph.linear, val %in% c(-3:4))
  graph.pop <- filter(pred.pop, list_class_ennahdha == "Ennahdha")
  

  # GRAPH BY POPULATION
  ggplot(data = graph.pop, aes(x = val, y = pred)) +
    # geom_linerange(aes(ymin = pred_l, ymax = pred_h, col = mun_type), size = .9, alpha = .4) +
    geom_ribbon(aes(ymin = pred_l, ymax = pred_h, fill = model), alpha = .05) +
    geom_line(aes(col = model), size = 3) +
    labs(x = "Relative strength in 2014 (standard deviations)",
         y = "Predicted\nprobability\nof FHL", color = "Log population", fill = "Log population",
         caption = "Predicted values based on a logit model controlling for party,\nwomen's turnout in 2014, urbanization rate, unemployment\nrate, and governorate fixed effects.")   +
    geom_hline(yintercept = .5, col = "black", size = .5, linetype = "dashed") +
    scale_y_continuous(limits = c(-1.2,1.2), labels = scales::percent_format(accuracy = 1)) +
    scale_x_continuous(breaks = seq(-3, 4, 1)) +
    scale_color_viridis(discrete = TRUE, end = .99) +
    scale_fill_viridis(discrete = TRUE, end = .99) +
    coord_cartesian(ylim=c(0, 1)) + lecs_theme
  ggsave("plot_pred_fhl_pop.pdf", plot = last_plot(),
         path = "figures", limitsize = F,
         width = 10, height = 6)
  
  ###### Check real population amounts
  real_pop <- list %>% dplyr::select(log_pop, sd_log_pop) %>%
    distinct() %>%
    mutate(pop = exp(log_pop))
  
    # +2 sd: total pop = 129693, Soukra, Ariana
    # +1 sd: total pop =  59578, Korba, Kalâa Kébira, Rades, Le Kef, Menzel Bourguiba
    # -1 sd: total pop =  10875, Lala, Hebira, Tibar
    # -2 sd: total pop =   4700, 	# El Masdour, Hazoua
  
  
#########################################################################   
#               ALTERNATIVE HYPOTHESIS                                  #   
#########################################################################       
  
  # Hypothesis: FHL LESS likely in parties/areas where voter bias against women is higher:
    # (a) municipalities where people are more religious (lower DI secular index)
    # (b) municipalities where people hold more patriarchal views about women in leadership (lower DI feminist index)
    # (c) for islamist parties (ennahdha)
  
  # Hypothesis: FHL MORE likely where there are more qualified women:
    # (d) municipalities that had more female representation in special delegations
    # (e) municipalities with more educated women
  
  ##### OLS MODELS [FOR APPENDIX TABLE A.3] ##### 
  
  # Create list of model names
  models.alt <- list()
  
  # Run models
  models.alt[[1]] <- lm(FHL ~ votes_per_nat_mean_2014_parl*list_class_ennahdha +
                            log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + gov_en, data = data.main)
 
  models.alt[[2]] <- lm(FHL ~ votes_per_nat_mean_2014_parl+ di_relig_pol_index +
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.main)
   
  models.alt[[3]] <- lm(FHL ~ votes_per_nat_mean_2014_parl+ di_women_pol_index +
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.main)
 
  models.alt[[4]] <- lm(FHL ~ votes_per_nat_mean_2014_parl+ sd_fem_per +
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.main)
  
  models.alt[[5]] <- lm(FHL ~ votes_per_nat_mean_2014_parl + edu_f_high + 
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.main)
  
  models.alt[[6]] <- lm(FHL ~ votes_per_nat_mean_2014_parl + edu_m_high + 
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.main)
  
  names(models.alt) <- c("Party Diff", "Secular Index", "Feminist Index", 
                         "Special Dels", "Women's Edu", "Men's Edu")
  
  # Cluster standard errors by municipality
  models.alt.cl <- list()
  models.alt.cl[[1]] <- cluster.se(models.alt[[1]], data.main$u_id)
  models.alt.cl[[2]] <- cluster.se(models.alt[[2]], data.main$u_id)
  models.alt.cl[[3]] <- cluster.se(models.alt[[3]], data.main$u_id)
  models.alt.cl[[4]] <- cluster.se(models.alt[[4]], data.main$u_id)
  models.alt.cl[[5]] <- cluster.se(models.alt[[5]], data.main$u_id)
  models.alt.cl[[6]] <- cluster.se(models.alt[[6]], data.main$u_id)

  # Calculate new F-stat with clustering
  check.f <- list()
  check.f[[1]] <- f.test(models.alt[[1]], data.main$u_id)
  check.f[[2]] <- f.test(models.alt[[2]], data.main$u_id)
  check.f[[3]] <- f.test(models.alt[[3]], data.main$u_id)
  check.f[[4]] <- f.test(models.alt[[4]], data.main$u_id)
  check.f[[5]] <- f.test(models.alt[[5]], data.main$u_id)
  check.f[[6]] <- f.test(models.alt[[6]], data.main$u_id)

  f <- NA
  p <- NA
  for(i in 1:length(check.f)){
    f[i] <- round(check.f[[i]]$F[2],3)
    p[i] <- round(check.f[[i]]$`Pr(>F)`[2],4)
  }
  p
  fstars <- data.frame(fval = f, pval = p, stars = rep(c("***"),length(models.alt))) %>%
    mutate(fstars = paste(fval, stars, sep = ""))
  
  # Output table
  stargazer(models.alt, 
            out = "tables/fhl_alt_ols.tex",
            title = "Alternative Explanations: Conservative Values and Candidate Supply",
            label = "fhl.alt.ols",
            table.placement = "H",
            type = "text", 
            se = models.alt.cl, 
            omit.stat = "f", 
            omit = c("gov_en", "mun_type"),
            font.size = "tiny",
            multicolumn = T,
            column.labels = names(models.alt),
            dep.var.labels = "Probability of an FHL",
            covariate.labels = c("Relative strength in 2014",
                                 "Secular Index (DI)",
                                 "Feminist Index (DI)",
                                 "Female SD members (%)",
                                 "Women w/ secondary edu or more (%)",
                                 "Men w/ secondary edu or more (%)",
                                 "Nidaa",
                                 "Third Party",
                                 "Log population",
                                 "Urbanization rate",
                                 "Unemployment rate",
                                 "Female turnout in 2014",
                                 "Relative strength x Nidaa",
                                 "Relative strength x Third Party",
                                 "Constant"),
            add.lines = list( #c("Controls", rep("Y", length(models.alt))),
                             c("Governorate FE", rep(c("Y"), length(models.alt))),
                             c("F Statistic", fstars$fstars)),
            notes.align = "l",
            notes = c("OLS models with standard errors clustered by municipality. All continuous variables (those except party) are standardized. Relative strength",
                      "is the party's voteshare at the municipal level in the 2014 parliamentary elections minus their average voteshare nationally. Both the secular",
                      "and feminist indices are municipal averages from a representative DI survey completed just after the 2018 municipal election and measure the'", 
                      "respondents' overall attitudes toward religion in politics (higher score = more secular) and women's leadership (higher score = more feminist).",
                      "Percent of female SD members is the percent of `Special Delegation` members--caretaker councils in place between the revolution and 2018 elections",
                      "---who were female.")
  ) 
  

  
  ##### LOGIT MODELS [FOR APPENDIX TABLE A.10] #####
  
  # Create list of model names
  models.alt.log <- list()
  
  # Run models
  models.alt.log[[1]] <- glm(FHL ~ votes_per_nat_mean_2014_parl*list_class_ennahdha +
                               log_pop  + pop_com_per + emp_unemploy_rate + 
                               turnout_per_pop_f_2014_parl + gov_en, data = data.main, family = "binomial")
  
  models.alt.log[[2]] <- glm(FHL ~ votes_per_nat_mean_2014_parl+ di_relig_pol_index +
                               log_pop  + pop_com_per + emp_unemploy_rate + 
                               turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                               gov_en, data = data.main, family = "binomial")
  
  models.alt.log[[3]] <- glm(FHL ~ votes_per_nat_mean_2014_parl+ di_women_pol_index +
                               log_pop  + pop_com_per + emp_unemploy_rate + 
                               turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                               gov_en, data = data.main, family = "binomial")
  
  models.alt.log[[4]] <- glm(FHL ~ votes_per_nat_mean_2014_parl+ sd_fem_per +
                               log_pop  + pop_com_per + emp_unemploy_rate + 
                               turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                               gov_en, data = data.main, family = "binomial")
  
  models.alt.log[[5]] <- glm(FHL ~ votes_per_nat_mean_2014_parl + edu_f_high + 
                               log_pop  + pop_com_per + emp_unemploy_rate + 
                               turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                               gov_en, data = data.main, family = "binomial")
  
  models.alt.log[[6]] <- glm(FHL ~ votes_per_nat_mean_2014_parl + edu_m_high + 
                               log_pop  + pop_com_per + emp_unemploy_rate + 
                               turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                               gov_en, data = data.main, family = "binomial")
  
  names(models.alt.log) <- c("Party Diff",  "Secular Index", "Feminist Index", 
                         "Special Dels", "Women's Edu", "Men's Edu")
  
  # Cluster standard errors by municipality
  models.alt.log.cl <- list()
  models.alt.log.cl[[1]] <- cluster.se(models.alt.log[[1]], data.main$u_id)
  models.alt.log.cl[[2]] <- cluster.se(models.alt.log[[2]], data.main$u_id)
  models.alt.log.cl[[3]] <- cluster.se(models.alt.log[[3]], data.main$u_id)
  models.alt.log.cl[[4]] <- cluster.se(models.alt.log[[4]], data.main$u_id)
  models.alt.log.cl[[5]] <- cluster.se(models.alt.log[[5]], data.main$u_id)
  models.alt.log.cl[[6]] <- cluster.se(models.alt.log[[6]], data.main$u_id)

  
#########################################################################   
#               APPENDIX: LOGIT MODELS                                  #   
#########################################################################   
  
  ###### FHL PLACEMENT: ALTERNATIVE EXPLANATIONS ######
  
  stargazer(models.alt.log, 
            out = "tables/fhl_alt_log.tex",
            title = "Logit Models: Alternative Explanations",
            label = "fhl.alt.log",
            table.placement = "H",
            type = "text", 
            se = models.alt.log.cl, 
            omit.stat = "f", 
            omit = c("gov_en", "mun_type"),
            font.size = "tiny",
            multicolumn = T,
            column.labels = names(models.alt.log),
            dep.var.labels = "Probability of an FHL",
            covariate.labels = c("Relative strength in 2014",
                                 "Secular Index (DI)",
                                 "Feminist Index (DI)",
                                 "Female SD members (%)",
                                 "Women w/ secondary edu or more (%)",
                                 "Men w/ secondary edu or more (%)",
                                 "Nidaa",
                                 "Third Party",
                                 "Log population",
                                 "Urbanization rate",
                                 "Unemployment rate",
                                 "Female turnout in 2014",
                                 "Relative strength x Nidaa",
                                 "Relative strength x Third Party",
                                 "Constant"),
            add.lines = list(#c("Controls", rep("Y", length(models.alt.log))),
                             c("Governorate FE", rep(c("Y"), length(models.alt.log))),
                             c("F Statistic", fstars$fstars)),
            notes.align = "l",
            notes = c("Logit models with standard errors clustered by municipality. All continuous variables (those except party) are standardized. Relative strength",
                      "is the party's voteshare at the municipal level in the 2014 parliamentary elections minus their average voteshare nationally. Both the secular",
                      "and feminist indices are municipal averages from a representative DI survey completed just after the 2018 municipal election and measure the'", 
                      "respondents' overall attitudes toward religion in politics (higher score = more secular) and women's leadership (higher score = more feminist).",
                      "Percent of female SD members is the percent of `Special Delegation` members--caretaker councils in place between the revolution and 2018 elections",
                      "---who were female.")
  ) 
  
#########################################################################   
#               APPENDIX: ALTERNATE STRENGTH MEASURES                   #   
#########################################################################    
  
  # Same results in terms of coefficients/directions, except squared
  # term on Quadratic non-FE model is significant at 10% level, but this 
  # disappears with FE.
  
  ###### OLS MODELS [FOR APPENDIX TABLE A.6] #####
  
  # Create list of model names
  strength.alt <- list()
  
  # Run models
  strength.alt[[1]] <-lm(FHL ~ votes_per_2014_parl + 
                           log_pop  + pop_com_per + emp_unemploy_rate + 
                           turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                           gov_en, data = data.main)
  
  strength.alt[[2]] <- lm(FHL ~ votes_per_2014_parl + votes_per_2014_parlsq + 
                            log_pop  + pop_com_per + emp_unemploy_rate + 
                            turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                            gov_en, data = data.main)
  
  strength.alt[[3]] <-lm(FHL ~ votes_per_dis_mean_2014_parl + 
                           log_pop  + pop_com_per + emp_unemploy_rate + 
                           turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                           gov_en, data = data.main)
  
  strength.alt[[4]] <- lm(FHL ~ votes_per_dis_mean_2014_parl + votes_per_dis_mean_2014_parlsq + 
                            log_pop  + pop_com_per + emp_unemploy_rate + 
                            turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                            gov_en, data = data.main)
  
  strength.alt[[5]] <-lm(FHL ~ list_margin3_2014_parl + 
                           log_pop  + pop_com_per + emp_unemploy_rate + 
                           turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                           gov_en, data = data.main)
  
  strength.alt[[6]] <- lm(FHL ~ list_margin3_2014_parl + list_margin3_2014_parlsq + 
                            log_pop  + pop_com_per + emp_unemploy_rate + 
                            turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                            gov_en, data = data.main)
  
  names(strength.alt) <- c("Voteshare", "Voteshare Quad", "District Mean", "District Mean Quad", "List Margin", "List Margin Quad")
  
  # Cluster standard errors by municipality
  strength.alt.cl <- list()
  strength.alt.cl[[1]] <- cluster.se(strength.alt[[1]], data.main$u_id)
  strength.alt.cl[[2]] <- cluster.se(strength.alt[[2]], data.main$u_id)
  strength.alt.cl[[3]] <- cluster.se(strength.alt[[3]], data.main$u_id)
  strength.alt.cl[[4]] <- cluster.se(strength.alt[[4]], data.main$u_id)
  strength.alt.cl[[5]] <- cluster.se(strength.alt[[5]], data.main$u_id)
  strength.alt.cl[[6]] <- cluster.se(strength.alt[[6]], data.main$u_id)
  
  # Calculate new F-stat with clustering
  check.f <- list()
  check.f[[1]] <- f.test(strength.alt[[1]], data.main$u_id)
  check.f[[2]] <- f.test(strength.alt[[2]], data.main$u_id)
  check.f[[3]] <- f.test(strength.alt[[3]], data.main$u_id)
  check.f[[4]] <- f.test(strength.alt[[4]], data.main$u_id)
  check.f[[5]] <- f.test(strength.alt[[5]], data.main$u_id)
  check.f[[6]] <- f.test(strength.alt[[6]], data.main$u_id)
  
  f <- NA
  p <- NA
  for(i in 1:length(check.f)){
    f[i] <- round(check.f[[i]]$F[2],3)
    p[i] <- round(check.f[[i]]$`Pr(>F)`[2],4)
  }
  p
  fstars <- data.frame(fval = f, pval = p, stars = rep(c("***"),length(strength.alt))) %>%
    mutate(fstars = paste(fval, stars, sep = ""))
  
  # Output table
  stargazer(strength.alt, 
            out = "tables/fhl_alt_strength_ols.tex",
            title = "Alternate Measures of Relative Strength",
            label = "strength.alt",
            type = "text", 
            column.labels = names(strength.alt),
            dep.var.labels = "Probability of an FHL",
            covariate.labels = c( "Voteshare in 2014",
                                  "Voteshare in 2014 sq",
                                  "Voteshare (dis. mean) in 2014",
                                  "Voteshare (dis. mean) in 2014 sq",
                                  "List margin of diff. in 2014",
                                  "List margin of diff. in 2014 sq",
                                  "Log population",
                                  "Constant"),
            se = strength.alt.cl, 
            omit.stat = "f", 
            omit = c("gov_en", "list_class_ennahdha", "mun_type","pop_com_per","emp_unemploy_rate","turnout_per_pop_f_2014_parl"),
            font.size = "scriptsize",
            add.lines = list(c("Controls", rep(c("Y"), length(strength.alt))),
                             c("Governorate FE", rep(c("Y"), length(strength.alt))),
                             c("F Statistic", fstars$fstars)),
            notes.align = "l",
            notes = c("OLS models with standard errors clustered by municipality, which include governorate fixed effects as well as female",
                      "turnout, log population, urbanization rate, unemployment rate, and party, all of which have effects and significance",
                      "levels comparable to the main models. Models 1-2 use the party's raw voteshare from the 2014 elections aggregated to",
                      "the municipality. Models 3-4 center this voteshare at the party's average voteshare across municipalities in the district.",
                      "Models 5-6 subtract the party's voteshare from the voteshare of the list that came in first place. Analysis excludes ",
                      "independents and parties that ran in 2018 but not 2014. Statistical significance and direction of coefficients is",
                      "consistent with logit models but OLS is shown here for ease of interpretation.")
  ) 
  

#########################################################################   
#               APPENDIX: ROBUSTNESS CHECKS                             #   
#########################################################################    
  
  ###### ALTERNATE SAMPLES AND SPECIFICATIONS 
  
  # Ennahdha and Nidaa only (no third parties)
  # Dropping biggest E/N strongholds (Tataouine, Monastir)
  # Drop largest cities (Tunis, Sfax, Sousse)
  # Using regional fixed effects
  
  ###### MODELS
  
  # Create list of model names
  robust.ols <- list()
  
  robust.ols[[1]] <- lm(FHL ~ votes_per_nat_mean_2014_parl + 
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.nn)
  
  robust.ols[[2]] <- lm(FHL ~ votes_per_nat_mean_2014_parl + 
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.nostrong)
  
  robust.ols[[3]] <- lm(FHL ~ votes_per_nat_mean_2014_parl + 
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          gov_en, data = data.nolarge)
  
  robust.ols[[4]] <- lm(FHL ~ votes_per_nat_mean_2014_parl + 
                          log_pop  + pop_com_per + emp_unemploy_rate + 
                          turnout_per_pop_f_2014_parl + list_class_ennahdha + 
                          reg, data = data.main)

  
  names(robust.ols) <- c("E/N Only", "Drop Strong", 
                         "Drop Large","Reg FE")
  
  
  # Cluster standard errors by municipality
  robust.ols.cl <- list()
  robust.ols.cl[[1]] <- cluster.se(robust.ols[[1]], data.nn$u_id)
  robust.ols.cl[[2]] <- cluster.se(robust.ols[[2]], data.nostrong$u_id)
  robust.ols.cl[[3]] <- cluster.se(robust.ols[[3]], data.nolarge$u_id)
  robust.ols.cl[[4]] <- cluster.se(robust.ols[[4]], data.main$u_id)

  # Calculate new F-stat with clustering
  check.f <- list()
  check.f[[1]] <- f.test(robust.ols[[1]], data.nn$u_id)
  check.f[[2]] <- f.test(robust.ols[[2]], data.nostrong$u_id)
  check.f[[3]] <- f.test(robust.ols[[3]], data.nolarge$u_id)
  check.f[[4]] <- f.test(robust.ols[[4]], data.main$u_id)

  f <- NA
  p <- NA
  for(i in 1:length(check.f)){
    f[i] <- round(check.f[[i]]$F[2],3)
    p[i] <- round(check.f[[i]]$`Pr(>F)`[2],4)
  }
  p
  fstars <- data.frame(fval = f, pval = p, stars = rep(c("***"),length(robust.ols))) %>%
    mutate(fstars = paste(fval, stars, sep = ""))
  
  # Output table
  stargazer(robust.ols, 
            out = "tables/fhl_robust_ols.tex",
            title = "Alternate Specifications of Main Models",
            label = "fhl.robust.ols",
            table.placement = "H",
            type = "text", 
            column.labels = names(robust.ols),
            dep.var.labels = "Probability of an FHL",
            covariate.labels = c("Relative strength in 2014",
                                 "Log population",
                                 "Urbanization rate",
                                 "Unemployment rate",
                                 "Female turnout in 2014",
                                 "Nidaa",
                                 "Third Party",
                                 "Constant"),
            se = robust.ols.cl, 
            omit.stat = "f", 
            omit = c("gov_en", "reg"),
            font.size = "tiny",
            add.lines = list(c("Governorate FE", rep("Y", 3), "N"),
                             c("Region FE", rep("N", 3), "Y"),
                             c("F Statistic", fstars$fstars)),
            notes.align = "l",
            notes = c("OLS models with standard errors clustered by municipality, and standardized continuous variables (relative strength, ", 
                      "log population, urbanization rate, unemployment rate, and female turnout). Model 1 drops Third Parties from the",
                      "analysis. Model 2 drops the two governorates where Ennahdha and Nidaa were the strongest (Tataouine and Monastir,",
                      "respectively). Model 3 drops the three largest municipalities in Tunisia (Tunis, Sfax, and Sousse). Model 4",
                      "uses regional fixed effects instead of governorate fixed effects.")
  ) 
  

  
  
